% This is file solves the models with the two different specifications 
% for the persistence of the productivity shock. It is being called by 
% 'FIG_12_Robustness_Theory_Persistence.m'

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Parameters and Steady-State Values %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


beta                    =   0.99;   % Discount factor.
chi                     =   0;      % Inverse of Frisch elasticity. 
sigma                   =   1;      % Inverse of EIS.

phi                     =   0.75;   % Prob. firms cannot adjust prices. 

%rho                     =   0.9;    % Persistence of natural rate shock. 

% Slope Phillips Curve:
gamma                   =   ((1-phi)*(1-phi*beta)/phi)*(sigma+chi);

T                       =   21+1;   % Number of periods.           

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             Creating IRFs                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

H                       =   3;      % Peg length. 

% Pre-allocating:
pi0                     =   zeros(T,1);
pi1                     =   zeros(T,1);
x0                      =   zeros(T,1);
x1                      =   zeros(T,1);
int                     =   zeros(T,1);
int1                    =   zeros(T,1);
a                       =   zeros(T,1);

% Initial values:
int(1:H)                =   -(1/beta-1)*ones(H,1);
int1(1:H)               =   -(1/beta-1)*ones(H,1);
x0(H)                   =   -(1/sigma)*int(H);
pi0(H)                  =   gamma*x0(H);
pi1(H)                  =   gamma*x1(H);

for jx = 1:H-1
    x0(H-jx)                =   x0(H-jx+1) - (1/sigma)*int(H-jx) + (1/sigma)*pi0(H-jx+1);
    pi0(H-jx)               =   gamma*x0(H-jx) + beta*pi0(H-jx+1);
end
y0                      =   x0;

% Now hit it with a one unit shock to a(t)
for jx = 1:T
    a(jx,1)                 =   rho^(jx-1)*1;
end
yf                      =   ((1+chi)/(sigma + chi))*a;
rf                      =   (sigma*(1+chi)*(rho-1)/(sigma + chi))*a;
x1(H)                   =   -(1/sigma)*int1(H) + (1/sigma)*rf(H);
int1(H+1:end)           =   rf(H+1:end);

% Creating IRFs:
for jx = 1:H-1
    x1(H-jx)                =   x1(H-jx+1) - (1/sigma)*int1(H-jx)...
                                + (1/sigma)*pi1(H-jx+1) + (1/sigma)*rf(H-jx);
    pi1(H-jx)               =   gamma*x1(H-jx) + beta*pi1(H-jx+1);
end

y1                      =   x1 + yf;
irfy3                   =   y1-y0;
irfpi3                  =   pi1-pi0;

%% Different Peg. 
H                       =   6;

% Pre-allocating:
pi0                     =   zeros(T,1);
pi1                     =   zeros(T,1);
x0                      =   zeros(T,1);
x1                      =   zeros(T,1);
int                     =   zeros(T,1);
int1                    =   zeros(T,1);
a                       =   zeros(T,1);

% Initial values:
int(1:H)                =   -(1/beta-1)*ones(H,1);
int1(1:H)               =   -(1/beta-1)*ones(H,1);
x0(H)                   =   -(1/sigma)*int(H);
pi0(H)                  =   gamma*x0(H);
pi1(H)                  =   gamma*x1(H);

for jx = 1:H-1
    x0(H-jx)                = x0(H-jx+1) - (1/sigma)*int(H-jx) + (1/sigma)*pi0(H-jx+1);
    pi0(H-jx)               = gamma*x0(H-jx) + beta*pi0(H-jx+1);
end
y0                      =   x0;

% Now hit it with a one unit shock to a(t)
for jx = 1:T
    a(jx,1)                 =   rho^(jx-1)*1;
end
yf                      =   ((1+chi)/(sigma + chi))*a;
rf                      =   (sigma*(1+chi)*(rho-1)/(sigma + chi))*a;
x1(H)                   =   -(1/sigma)*int1(H) + (1/sigma)*rf(H);
int1(H+1:end)           =   rf(H+1:end);

% Creating IRFs:
for jx = 1:H-1
    x1(H-jx)                =   x1(H-jx+1) - (1/sigma)*int1(H-jx)...
                                + (1/sigma)*pi1(H-jx+1) + (1/sigma)*rf(H-jx);
    pi1(H-jx)               =   gamma*x1(H-jx) + beta*pi1(H-jx+1);
end


y1                      =   x1 + yf;
irfy6                   =   y1-y0;
irfpi6                  =   pi1-pi0;

%% Different Peg. 
H                       =   10;

% Pre-allocating:
pi0                     =   zeros(T,1);
pi1                     =   zeros(T,1);
x0                      =   zeros(T,1);
x1                      =   zeros(T,1);
int                     =   zeros(T,1);
int1                    =   zeros(T,1);
a                       =   zeros(T,1);

% Initial values:
int(1:H)                =   -(1/beta-1)*ones(H,1);
int1(1:H)               =   -(1/beta-1)*ones(H,1);
x0(H)                   =   -(1/sigma)*int(H);
pi0(H)                  =   gamma*x0(H);
pi1(H)                  =   gamma*x1(H);

for jx = 1:H-1
    x0(H-jx)                =   x0(H-jx+1) - (1/sigma)*int(H-jx) + (1/sigma)*pi0(H-jx+1);
    pi0(H-jx)               = gamma*x0(H-jx) + beta*pi0(H-jx+1);
end
y0                      =   x0;

% Now hit it with a one unit shock to a(t)
for jx = 1:T
    a(jx,1)                 =   rho^(jx-1)*1;
end
yf                      =   ((1+chi)/(sigma + chi))*a;
rf                      =   (sigma*(1+chi)*(rho-1)/(sigma + chi))*a;
x1(H)                   =   -(1/sigma)*int1(H) + (1/sigma)*rf(H);
int1(H+1:end)           =   rf(H+1:end);

% Creating IRFs:
for jx = 1:H-1
    x1(H-jx)                =   x1(H-jx+1) - (1/sigma)*int1(H-jx)...
                                + (1/sigma)*pi1(H-jx+1) + (1/sigma)*rf(H-jx);
    pi1(H-jx)               =   gamma*x1(H-jx) + beta*pi1(H-jx+1);
end


y1                      =   x1 + yf;
irfy10                  =   y1-y0;
irfpi10                 =   pi1-pi0;

pif                     =   zeros(T,1);

